function [data,data_trans] = generate_data_simult(N,J,mu_x,sigma_x,mu_z,sigma_z,alpha,beta,...
                                           mu_c,sigma_c)

% draw data

x = randn(N,J)*sigma_x + mu_x;
z = randn(N,J)*sigma_z + mu_z;

c = randn(N,1)*sigma_c + mu_c;
epsilon = -evrnd(0,1,[N,J]);

choice = zeros(N,1);
open = zeros(N,J);
VU = zeros(N,J);
U = zeros(N,J);

for i=1:N
       
    c_temp = c(i);
    
    VU(i,:) = alpha*x(i,:) + epsilon(i,:);
    
    [~,ind_VU] = sort(VU(i,:),'descend');
    
    open(i,:) = (VU(i,:)>c_temp);
    open(i,ind_VU(1)) = 1;       % good with highest VU is always searched
    
    U(i,:) = alpha*x(i,:) + beta*z(i,:) + epsilon(i,:);
    
    [~,max_temp] = max(U(i,logical(open(i,:))));
    ind_temp = find(open(i,:));
    choice(i) = ind_temp(max_temp);
            
end

data.x = x;
data.z = z;
data.c = c;
data.epsilon = epsilon;
data.choice = choice;
data.open = open;
data.VU = VU;
data.U = U;
% data.y = zeros(size(choice));
% data.y(choice==1) = 1;

data = struct2dataset(data);

data_trans = transform_data_01(data,J);

end

